Wunsch I · Perfekte Daten
In der Vorlesung 🧑🏫 am 19. November 2020 gab es eine Fee 🧚 die Wünsche 🧞 erfüllt.

Den Wunsch 🪄 nach perfekten Daten 📈 erfüllt sie nicht 😥.
Die Fee 🧚 überlegt, nutzt R 👩💻 und schafft fast perfekte Daten 📈.
param_const <- 500
param_alter <- 25
n <- 1000
sd_einkommen <- 100
dt_raw <-
tibble(
id = 1:n,
alter = seq(20, 60, length.out = n) %>% round(2),
einkommen_linear = round(param_const + param_alter*alter)
) %>%
mutate(einkommen = round(einkommen_linear + rnorm(n, sd = sd_einkommen)))
\[ Einkommen_{ID} = 500 + 25 \cdot Alter_{ID} \: + ( \: Fehler_{ID} \: )\]
ggplot(dt_raw, aes(x = alter, y = einkommen)) +
geom_point(alpha = 0.2) +
geom_smooth(method = lm) +
theme_stata()

Wunsch II · Diagnose für SoWis
Die Studierenden 🙇🙇🙇🙇 hatten einen eigenen Wunsch 🪄 an die Fee 🧚.
Kannst Du uns Regressions-Diagnose 📈 mit einer sozialwissenschaftlichen Geschichte 📚 erklären?

Daten Fee
Wir stellen uns Postbeamte 📯 vor die regelmäßig Lohnerhöhungen 💰 strikt nach dem Senioritätsprinzip 🥸 erhalten.
Wir haben also ein Modell 📈 bei dem das Einkommen 💰 vom Alter 🧑👴 abhängt.

Das wahre Einkommen 💰 (einkommen_linear) ist uns bekannt und unsere beobachteten Daten 🕵️ (einkommen) weichen davon normal-verteilt 🧮 ab
So hat die Fee 🧚 uns die Daten 📊 gegeben.
ggplot(dt_raw, aes(x = alter, y = einkommen)) +
geom_point(alpha = 0.2) +
geom_jitter(aes(y = einkommen_linear),
width = 0.1, alpha = 0.5,
colour = "darkgreen", shape = ".") +
theme_stata()

DT::datatable(dt_raw)
Modell (fast) perfekt
Wir bekommen schöne Ergebnisse bei der Schätzung des Modells. 😄
fit <- lm(einkommen ~ alter, data = dt_raw)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Annahmen (1-7)
- korrekte Spezifikation (Linearität)
- Erwartungswert der Residuen (lokal) = 0
- Normalverteilte Residuen
- Varianzhomogenität (gleichmäßige Streuung Residuen)
- Berücksichtigung einflussreicher Fälle (keine Ausreisser)
- keine Multikollinearität (geringe Korrelation unabhängiger Variablen)
- keine Autokorrelation der Residuen (insbes. Zeitreihen)
Korrekte Spezifikation
Dritt-Variablen

Briefträger ✉️ arbeiten oft nur halbtags 🕐 und verdienen 💰 dann auch nur die Hälfte.
dt <-
dt_raw %>%
mutate(einkommen = if_else(id %% 3 == 0, einkommen/2 + rnorm(n, sd = 50), einkommen),
halbtags = if_else(id %% 3 == 0, "ja", "nein") %>% fct_relevel("ja", after = Inf))
ggplot(dt, aes(x = alter, y = einkommen, colour = halbtags)) +
geom_point(alpha = 0.2) +
theme_stata()

Ohne Dritt-Variable 😟
fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit, bw=50)

Mit Dritt-Variable 🙂
fit <- lm(einkommen ~ alter + halbtags, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Noch besser mit Dritt-Variable und Interaktions-Effekt 😃
fit <- lm(einkommen ~ alter*halbtags, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Varianzhomogenität

Bei der Post 📯 haben hohe Beamte und Minster ein höheres Einkommen 💰 als einfache Beamte. Aber nicht alle werden hohe Beamte 🧑💼 oder Post-Minister 🏤.
dt <- dt_raw %>% mutate(einkommen = einkommen_linear + id/500 * rnorm(n, sd = sd_einkommen))
ggplot(dt, aes(x = alter, y = einkommen)) +
geom_point(alpha = 0.2) +
geom_jitter(aes(y = einkommen_linear),
width = 0.1, alpha = 0.5,
colour = "darkgreen", shape = ".") +
theme_stata()

Im Standard-Modell erhalten wir verzerrte Standard-Fehler. 😟
fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Besser ist ein Modell mit robusten 💪 Standard-Fehlern. 😄
fit <- lm_robust(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance(fit) %>% mutate_if(is.numeric, round, 2)
Multikollinearität

Ältere 👴👵 Postbeamte 📯 sind kleiner als Jüngere 👧👦.
dt <- dt_raw %>% mutate(groesse_cm = round(seq(180, 160, length.out = n) + rnorm(n, sd = 2.5), 3))
ggplot(dt, aes(x = alter, y = groesse_cm)) +
geom_point(alpha = 0.2) +
geom_smooth(method = lm) +
theme_stata()

Es gibt dann auch einen Zusammenhang 📉 zwischen groesse_cm 👵 und einkommen 💰, natürlich keinen kausalen ⚙️.
ggplot(dt, aes(x = groesse_cm, y = einkommen)) +
geom_point(alpha = 0.2) +
geom_smooth(method = lm) +
theme_stata()

fit <- lm(einkommen ~ groesse_cm, data = dt)
tidy_lm(fit)
glance_lm(fit)
Nehmen wir alter 👵 und groesse_cm 📐 in einem Modell 📉 auf haben wir Multikollinearität 🔗 😟
fit <- lm(einkommen ~ alter + groesse_cm, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Multikollinearität im Modell 📉 erkennen 🕵️ wir auch an der sehr starken Korrellation der unabhängigen Variablen.
dt %>%
select(alter, einkommen) %>%
cor() %>%
as_tibble() %>%
mutate_all(round, 2)
Ausreisser

Bei der Datenübertragung 📡 gab es einen Fehler ⚠️ und die letzte Ziffer 🔢 wurde bei einigen Gehaltsangaben 💰 entfernt.
id_random <- runif(80, min = 1, max = n) %>% round()
dt <-
dt_raw %>%
mutate(einkommen = if_else(id %in% id_random, round(einkommen/10), einkommen))
ggplot(dt, aes(x = alter, y = einkommen)) +
geom_point(alpha = 0.2) +
theme_stata()

Wir erhalten eine verzerrte Schätzung. 😟
fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit, bw=50)

Die falschen Beobachtungen müssen wir entfernen (oder korrigieren). 😃
dt_tmp <- dt %>% filter(einkommen > 500)
fit <- lm(einkommen ~ alter, data = dt_tmp)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Autokorrelation

Die Postbeamten 📯 wurden 2010 und 2020 📅 nach ihrem Einkommen 💰 befragt. ☎️
dt_tmp <-
dt_raw %>%
mutate(einkommen = round(einkommen - param_alter*10 + rnorm(n, sd = 10)),
alter = alter - 10,
befragung = 2010) %>%
filter(alter >= 20)
dt <-
dt_raw %>%
mutate(befragung = 2020) %>%
rbind(dt_tmp) %>%
arrange(id, befragung) %>%
relocate(befragung, .after = id)
DT::datatable(dt %>% filter(id > 500)) # Daten ab ID 500 um Panel-Struktur zu verdeutlichen
Postbeamte 📯 die 2010 📅 älter als 50 waren sind inzwischen in Pension 👵 und haben an der 2020 Befragung ☎️ nicht mehr teilgenommen.
dt_tmp <- dt %>% mutate(befragung = factor(befragung))
ggplot(dt_tmp, aes(x = alter, y = einkommen, colour = befragung)) +
geom_point(alpha = 0.2) +
theme_stata()

Im Modell (n = 1750) sind Beobachtungen (und Residuen) nicht mehr unabhängig. 😟 Wir haben ja viele zweimal befragt. ☎️
fit <- lm(einkommen ~ alter, data = dt_tmp)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

LS0tCnRpdGxlOiAiUmVncmVzc2lvbnMtRGlhZ25vc2Ugwrcg8J+UjSIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKICAgIHRvYzogeWVzCiAgICB0b2NfZmxvYXQ6IHllcwotLS0KCgpgYGB7ciwgbWVzc2FnZT1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCgpsaWJyYXJ5KGJyb29tKSAgICAgIyB0aWR5IG1vZGVsIGRhdGEKbGlicmFyeShwYXRjaHdvcmspICMgY29tYmluZSBnZ3Bsb3RzCmxpYnJhcnkoZ2d0aGVtZXMpICAjIGdncGxvdCBTdGF0YSB0aGVtZQpsaWJyYXJ5KGVzdGltYXRyKSAgIyByb2J1c3Qgc3RhbmRhcmQtZXJyb3JzCmBgYAoKCmBgYHtyLCBtZXNzYWdlPUZBTFNFfQojIyAgdGlkeSBsbSBtb2RlbCByZXN1bHRzIHdpdGggYnJvb20gKGFuZCBnZ3Bsb3QpCgp0aWR5X2xtIDwtIGZ1bmN0aW9uKGZpdCkgewogIHRpZHkoZml0KSAlPiUgbXV0YXRlX2lmKGlzLm51bWVyaWMsIHJvdW5kLCAyKQp9CgpnbGFuY2VfbG0gPC0gZnVuY3Rpb24oZml0KSB7CiAgZ2xhbmNlKGZpdCkgJT4lIAogICAgc2VsZWN0KHIuc3F1YXJlZCwgYWRqLnIuc3F1YXJlZCwgc3RhdGlzdGljLCBwLnZhbHVlLCBkZiwgbm9icykgJT4lIAogICAgbXV0YXRlX2lmKGlzLm51bWVyaWMsIHJvdW5kLCAyKQp9CgpwbG90X2xtIDwtIGZ1bmN0aW9uKGZpdCwgICBidyA9IDI1KSB7CiAgZHRfcGw8LSAgCiAgICBhdWdtZW50KGZpdCkgJT4lIAogICAgbXV0YXRlKGVpbmtvbW1lbl9nZXNjaGFldHp0ID0gLmZpdHRlZCwKICAgICAgICAgICByZXNpZHVlbiA9IC5yZXNpZCkKCiAgcGwxIDwtIAogICAgZ2dwbG90KGR0X3BsLCBhZXMoeCA9IHJlc2lkdWVuKSkgKwogICAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSBidywgY29sb3VyID0gImJsYWNrIiwgZmlsbCA9ICJncmV5IikgKwogICAgc3RhdF9mdW5jdGlvbigKICAgICAgZnVuID0gZnVuY3Rpb24oeCkgZG5vcm0oeCwgbWVhbiA9IG1lYW4oZHRfcGwkcmVzaWR1ZW4pLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzZCA9IHNkKGR0X3BsJHJlc2lkdWVuKSkKICAgICAgICAgICAgICAgICAgICAgICAgKiBidyAqIG5yb3coZHRfcGwpLAogICAgICBjb2xvdXIgPSAiYmx1ZSIKICAgICkgKwogICB0aGVtZV9zdGF0YSgpCgogIHBsMiA8LSAKICAgIGdncGxvdChkYXRhID0gZHRfcGwsIGFlcyh4ID0gZWlua29tbWVuX2dlc2NoYWV0enQsIHkgPSByZXNpZHVlbikpICsKICAgIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICAgIGdlb21fc21vb3RoKHNlID0gRkFMU0UpICsKICAgIHRoZW1lX3N0YXRhKCkKCiAgcGwxICsgcGwyICsgcGxvdF9sYXlvdXQod2lkdGhzID0gYygyLCAzKSkKfQpgYGAKCgojIFd1bnNjaCBJIMK3IFBlcmZla3RlIERhdGVuCgpJbiBkZXIgVm9ybGVzdW5nIPCfp5HigI3wn4+rIGFtIDE5LiBOb3ZlbWJlciAyMDIwIGdhYiBlcyBlaW5lIEZlZSDwn6eaIGRpZSBXw7xuc2NoZSDwn6eeIGVyZsO8bGx0LgoKIVtdKGJpbGRlci9mYWlyLmpwZykKCkRlbiBXdW5zY2gg8J+qhCBuYWNoIHBlcmZla3RlbiBEYXRlbiDwn5OIIGVyZsO8bGx0IHNpZSBuaWNodCDwn5ilLgogCkRpZSBGZWUg8J+nmiDDvGJlcmxlZ3QsIG51dHp0IFtfX1JfX10oaHR0cHM6Ly9lZHVjYXRpb24ucnN0dWRpby5jb20vbGVhcm4vYmVnaW5uZXIvKSDwn5Gp4oCN8J+SuyB1bmQgc2NoYWZmdCBmYXN0IHBlcmZla3RlIERhdGVuIPCfk4guCgpgYGB7cn0KcGFyYW1fY29uc3QgPC0gNTAwCnBhcmFtX2FsdGVyIDwtIDI1CgpuIDwtIDEwMDAKc2RfZWlua29tbWVuIDwtICAxMDAKCmR0X3JhdyA8LSAKICB0aWJibGUoIAogICAgaWQgPSAxOm4sCiAgICBhbHRlciA9IHNlcSgyMCwgNjAsIGxlbmd0aC5vdXQgPSBuKSAlPiUgcm91bmQoMiksCiAgICBlaW5rb21tZW5fbGluZWFyID0gcm91bmQocGFyYW1fY29uc3QgKyBwYXJhbV9hbHRlciphbHRlcikKICAgICkgJT4lIAogIG11dGF0ZShlaW5rb21tZW4gPSByb3VuZChlaW5rb21tZW5fbGluZWFyICsgcm5vcm0obiwgc2QgPSBzZF9laW5rb21tZW4pKSkKYGBgCgoKJCQgRWlua29tbWVuX3tJRH0gPSA1MDAgKyAyNSBcY2RvdCBBbHRlcl97SUR9IFw6ICsgKCBcOiBGZWhsZXJfe0lEfSBcOiApJCQKCmBgYHtyfQpnZ3Bsb3QoZHRfcmF3LCBhZXMoeCA9IGFsdGVyLCB5ID0gZWlua29tbWVuKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSBsbSkgKwogIHRoZW1lX3N0YXRhKCkKYGBgCgoKIyBXdW5zY2ggSUkgwrcgRGlhZ25vc2UgZsO8ciBTb1dpcwoKRGllIFN0dWRpZXJlbmRlbiDwn5mH8J+Zh/CfmYfwn5mHIGhhdHRlbiBlaW5lbiBlaWdlbmVuIFd1bnNjaCDwn6qEIGFuIGRpZSBGZWUg8J+nmi4KCkthbm5zdCBEdSB1bnMgUmVncmVzc2lvbnMtRGlhZ25vc2Ug8J+TiCBtaXQgZWluZXIgc296aWFsd2lzc2Vuc2NoYWZ0bGljaGVuIEdlc2NoaWNodGUg8J+TmiBlcmtsw6RyZW4/CgohW10oYmlsZGVyL2ZhaXIuanBnKQoKIyMgRGF0ZW4gRmVlCgpXaXIgc3RlbGxlbiB1bnMgUG9zdGJlYW10ZSDwn5OvIHZvciBkaWUgcmVnZWxtw6TDn2lnIExvaG5lcmjDtmh1bmdlbiDwn5KwIHN0cmlrdCBuYWNoIGRlbSBTZW5pb3JpdMOkdHNwcmluemlwIPCfpbggZXJoYWx0ZW4uIAoKV2lyIGhhYmVuIGFsc28gZWluIE1vZGVsbCDwn5OIIGJlaSBkZW0gZGFzIF9fRWlua29tbWVuX18g8J+SsCB2b20gX19BbHRlcl9fIPCfp5Hwn5G0IGFiaMOkbmd0LgoKIVtdKGJpbGRlci9wb3N0YmVhbXRlLmpwZykKCkRhcyB3YWhyZSBFaW5rb21tZW4g8J+SsCAoYGVpbmtvbW1lbl9saW5lYXJgKSBpc3QgdW5zIGJla2FubnQgdW5kIHVuc2VyZSBiZW9iYWNodGV0ZW4gRGF0ZW4g8J+Vte+4jyAoYGVpbmtvbW1lbmApIHdlaWNoZW4gZGF2b24gbm9ybWFsLXZlcnRlaWx0IPCfp64gYWIKClNvIGhhdCBkaWUgRmVlIPCfp5ogdW5zIGRpZSBEYXRlbiDwn5OKIGdlZ2ViZW4uCgpgYGB7cn0KZ2dwbG90KGR0X3JhdywgYWVzKHggPSBhbHRlciwgeSA9IGVpbmtvbW1lbikpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC4yKSArCiAgZ2VvbV9qaXR0ZXIoYWVzKHkgPSBlaW5rb21tZW5fbGluZWFyKSwKICAgICAgICAgICAgICB3aWR0aCA9IDAuMSwgYWxwaGEgPSAwLjUsCiAgICAgICAgICAgICAgY29sb3VyID0gImRhcmtncmVlbiIsIHNoYXBlID0gIi4iKSArCiAgdGhlbWVfc3RhdGEoKQoKRFQ6OmRhdGF0YWJsZShkdF9yYXcpCmBgYAoKIyMgTW9kZWxsIChmYXN0KSBwZXJmZWt0CgpXaXIgYmVrb21tZW4gc2Now7ZuZSBFcmdlYm5pc3NlIGJlaSBkZXIgU2Now6R0enVuZyBkZXMgTW9kZWxscy4g8J+YhAoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0X3JhdykKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCgojIyBBbm5haG1lbiAoMS03KQoKMS4ga29ycmVrdGUgU3BlemlmaWthdGlvbiAoTGluZWFyaXTDpHQpCjIuIEVyd2FydHVuZ3N3ZXJ0IGRlciBSZXNpZHVlbiAobG9rYWwpID0gMAozLiBOb3JtYWx2ZXJ0ZWlsdGUgUmVzaWR1ZW4KNC4gVmFyaWFuemhvbW9nZW5pdMOkdCAoZ2xlaWNobcOkw59pZ2UgU3RyZXV1bmcgUmVzaWR1ZW4pCjUuIEJlcsO8Y2tzaWNodGlndW5nIGVpbmZsdXNzcmVpY2hlciBGw6RsbGUgKGtlaW5lIEF1c3JlaXNzZXIpCjYuIF9rZWluZV8gTXVsdGlrb2xsaW5lYXJpdMOkdCAoZ2VyaW5nZSBLb3JyZWxhdGlvbiB1bmFiaMOkbmdpZ2VyIFZhcmlhYmxlbikKNy4gX2tlaW5lXyBBdXRva29ycmVsYXRpb24gZGVyIFJlc2lkdWVuIChpbnNiZXMuIFplaXRyZWloZW4pCgoKIyMgS29ycmVrdGUgU3BlemlmaWthdGlvbgoKIyMjIERyaXR0LVZhcmlhYmxlbgoKIVtdKGJpbGRlci9icmllZnRyYWVnZXItZnJhdWVuLmpwZykKCkJyaWVmdHLDpGdlciDinInvuI8gYXJiZWl0ZW4gb2Z0IG51ciBoYWxidGFncyDwn5WQIHVuZCB2ZXJkaWVuZW4g8J+SsCBkYW5uIGF1Y2ggbnVyIGRpZSBIw6RsZnRlLgoKYGBge3J9CmR0IDwtCiAgZHRfcmF3ICU+JSAKICBtdXRhdGUoZWlua29tbWVuID0gaWZfZWxzZShpZCAlJSAzID09IDAsIGVpbmtvbW1lbi8yICsgcm5vcm0obiwgc2QgPSA1MCksIGVpbmtvbW1lbiksCiAgICAgICAgIGhhbGJ0YWdzID0gaWZfZWxzZShpZCAlJSAzID09IDAsICJqYSIsICJuZWluIikgJT4lIGZjdF9yZWxldmVsKCJqYSIsIGFmdGVyID0gSW5mKSkKCmdncGxvdChkdCwgYWVzKHggPSBhbHRlciwgeSA9IGVpbmtvbW1lbiwgY29sb3VyID0gaGFsYnRhZ3MpKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuMikgKwogIHRoZW1lX3N0YXRhKCkKYGBgCgpPaG5lIERyaXR0LVZhcmlhYmxlIPCfmJ8KCmBgYHtyfQpmaXQgPC0gbG0oZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0LCBidz01MCkKYGBgCgpNaXQgRHJpdHQtVmFyaWFibGUg8J+ZggoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciArIGhhbGJ0YWdzLCBkYXRhID0gZHQpCnRpZHlfbG0oZml0KQpnbGFuY2VfbG0oZml0KQpwbG90X2xtKGZpdCkKYGBgCgpOb2NoIGJlc3NlciBtaXQgRHJpdHQtVmFyaWFibGUgdW5kIEludGVyYWt0aW9ucy1FZmZla3Qg8J+YgwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlcipoYWxidGFncywgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKcGxvdF9sbShmaXQpCmBgYAoKIyMjIEZ1bmt0aW9uYWxlIEZvcm0KCiFbXShiaWxkZXIvYnVuZGVzcG9zdC5qcGcpCgpCZWkgZGVyIEJ1bmRlc3Bvc3Qg8J+TryB3dXJkZSBlaW5lIEJlemFobHVuZyDwn5KwIG5hY2ggSHVtYW5rYXBpdGFsIGVpbmdlZsO8aHJ0LiAKCkVzIHd1cmRlIGVya2FubnQsIGRhc3MgYmVpIGp1bmdlbiBCZWFtdGVuIPCfp5EgZGllIFByb2R1a3Rpdml0w6R0IHN0w6Rya2VyIHfDpGNoc3QgYWxzIGJlaSDDhGx0ZXJlbiDwn6eTLiBEYWhlciBlcmhhbHRlbiBzaWUgaMO2aGVyZSBMb2huenV3w6RjaHNlIPCfkrAuCgpgYGB7cn0KZHQgPC0KICBkdF9yYXcgJT4lIAogIG11dGF0ZShlaW5rb21tZW5fcG9seW5vbSA9IC0zMDAgKyA4MCphbHRlciAtIDAuNzI1KmFsdGVyXjIsCiAgICAgICAgIGVpbmtvbW1lbiA9IGVpbmtvbW1lbl9wb2x5bm9tICsgcm5vcm0obiwgc2QgPSBzZF9laW5rb21tZW4pKQoKZ2dwbG90KGR0LCBhZXMoeCA9IGFsdGVyLCB5ID0gZWlua29tbWVuKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICAgIGdlb21faml0dGVyKGFlcyh5ID0gZWlua29tbWVuX3BvbHlub20pLAogICAgICAgICAgICAgICAgd2lkdGggPSAwLjEsIGFscGhhID0gMC41LAogICAgICAgICAgICAgICAgY29sb3VyID0gImRhcmtibHVlIiwgc2hhcGUgPSAiLiIpICsKICB0aGVtZV9zdGF0YSgpCmBgYApFaW4gbGluZWFyZXMgTW9kZWxsIGlzdCBmYWxzY2ggc3BlemlmaXppZXJ0LiDwn5ifIERpZSBSZWdyZXNzaW9ucy1EaWFnbm9zZSB6ZWlndCBkaWVzLiDimJ3vuI8KCmBgYHtyfQpmaXQgPC0gbG0oZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCktvcnJla3QgaXN0IGVpbiBNb2RlbGwgbWl0IGVpbmVtIFBvbHlub20gKGhpZXIgZWluZW0gcXVhZHJhdGlzY2hlbSBUZXJtKS4g8J+YgwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciArIEkoYWx0ZXJeMiksIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCiggRGVuIFIgSGlud2VpcyBhdWYgTXVsdGktS29sbGluZWFyaXTDpHQgaWdub3JpZXJlbiB3aXIgZWluZmFjaC4gTmF0w7xybGljaCBnaWJ0IGVzIGVpbmVuIFp1c2FtbWVuaGFuZyB6d2lzY2hlbiBgYWx0ZXJgIHVuZCBkZW0gcXVhZHJpZXJ0ZW4gYGFsdGVyXjJgLiApCgoKIyMgVmFyaWFuemhvbW9nZW5pdMOkdAoKIVtdKGJpbGRlci9taW5pc3Rlci5qcGcpCgpCZWkgZGVyIFBvc3Qg8J+TryBoYWJlbiBob2hlIEJlYW10ZSB1bmQgTWluc3RlciBlaW4gaMO2aGVyZXMgRWlua29tbWVuIPCfkrAgYWxzIGVpbmZhY2hlIEJlYW10ZS4gQWJlciBuaWNodCBhbGxlIHdlcmRlbiBob2hlIEJlYW10ZSDwn6eR4oCN8J+SvCBvZGVyIFtQb3N0LU1pbmlzdGVyXShodHRwczovL2RlLndpa2lwZWRpYS5vcmcvd2lraS9DaHJpc3RpYW5fU2Nod2Fyei1TY2hpbGxpbmcpIPCfj6QuCgpgYGB7cn0KZHQgPC0gZHRfcmF3ICU+JSBtdXRhdGUoZWlua29tbWVuID0gZWlua29tbWVuX2xpbmVhciArIGlkLzUwMCAqIHJub3JtKG4sIHNkID0gc2RfZWlua29tbWVuKSkKCmdncGxvdChkdCwgYWVzKHggPSBhbHRlciwgeSA9IGVpbmtvbW1lbikpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC4yKSArCiAgICBnZW9tX2ppdHRlcihhZXMoeSA9IGVpbmtvbW1lbl9saW5lYXIpLAogICAgICAgICAgICAgICAgd2lkdGggPSAwLjEsIGFscGhhID0gMC41LAogICAgICAgICAgICAgICAgY29sb3VyID0gImRhcmtncmVlbiIsIHNoYXBlID0gIi4iKSArCiAgdGhlbWVfc3RhdGEoKQpgYGAKCkltIFN0YW5kYXJkLU1vZGVsbCBlcmhhbHRlbiB3aXIgdmVyemVycnRlIFN0YW5kYXJkLUZlaGxlci4g8J+YnwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKcGxvdF9sbShmaXQpCmBgYAoKQmVzc2VyIGlzdCBlaW4gTW9kZWxsIG1pdCBfcm9idXN0ZW5fIPCfkqogU3RhbmRhcmQtRmVobGVybi4g8J+YhAoKYGBge3J9CmZpdCA8LSBsbV9yb2J1c3QoZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZShmaXQpICU+JSBtdXRhdGVfaWYoaXMubnVtZXJpYywgcm91bmQsIDIpIApgYGAKCiMjIE11bHRpa29sbGluZWFyaXTDpHQKCiFbXShiaWxkZXIvYnJpZWZ0cmFlZ2VyLWdydXBwZS5qcGcpCgrDhGx0ZXJlIPCfkbTwn5G1IFBvc3RiZWFtdGUg8J+TryBzaW5kIGtsZWluZXIgYWxzIErDvG5nZXJlIPCfkafwn5GmLgoKYGBge3J9CmR0IDwtIGR0X3JhdyAlPiUgbXV0YXRlKGdyb2Vzc2VfY20gPSByb3VuZChzZXEoMTgwLCAxNjAsIGxlbmd0aC5vdXQgPSBuKSArIHJub3JtKG4sIHNkID0gMi41KSwgMykpCgpnZ3Bsb3QoZHQsIGFlcyh4ID0gYWx0ZXIsIHkgPSBncm9lc3NlX2NtKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSBsbSkgKwogIHRoZW1lX3N0YXRhKCkKYGBgCkVzIGdpYnQgZGFubiBhdWNoIGVpbmVuIFp1c2FtbWVuaGFuZyDwn5OJIHp3aXNjaGVuIGBncm9lc3NlX2NtYCDwn5G1IHVuZCBgZWlua29tbWVuYCDwn5KwLCBuYXTDvHJsaWNoIGtlaW5lbiBrYXVzYWxlbiDimpnvuI8uCgpgYGB7cn0KZ2dwbG90KGR0LCBhZXMoeCA9IGdyb2Vzc2VfY20sIHkgPSBlaW5rb21tZW4pKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuMikgKwogIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtKSArCiAgdGhlbWVfc3RhdGEoKQpgYGAKCmBgYHtyfQpmaXQgPC0gbG0oZWlua29tbWVuIH4gZ3JvZXNzZV9jbSwgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKYGBgCgpOZWhtZW4gd2lyIGBhbHRlcmAg8J+RtSB1bmQgYGdyb2Vzc2VfY21gIPCfk5AgaW4gZWluZW0gTW9kZWxsIPCfk4kgYXVmIGhhYmVuIHdpciBNdWx0aWtvbGxpbmVhcml0w6R0IPCflJcg8J+YnwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciArIGdyb2Vzc2VfY20sIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCk11bHRpa29sbGluZWFyaXTDpHQgaW0gTW9kZWxsIPCfk4kgIGVya2VubmVuIPCflbXvuI8gd2lyIGF1Y2ggIGFuIGRlciBzZWhyIHN0YXJrZW4gS29ycmVsbGF0aW9uIGRlciB1bmFiaMOkbmdpZ2VuIFZhcmlhYmxlbi4KCmBgYHtyfQpkdCAlPiUgCiAgc2VsZWN0KGFsdGVyLCBlaW5rb21tZW4pICU+JSAKICBjb3IoKSAlPiUgCiAgYXNfdGliYmxlKCkgJT4lIAogIG11dGF0ZV9hbGwocm91bmQsIDIpCmBgYAoKCiMjIEF1c3JlaXNzZXIKCiFbXShiaWxkZXIvZmVybm1lbGRlYW10LmpwZykKCkJlaSBkZXIgRGF0ZW7DvGJlcnRyYWd1bmcg8J+ToSBnYWIgZXMgZWluZW4gRmVobGVyIOKaoO+4jyB1bmQgZGllIGxldHp0ZSBaaWZmZXIg8J+UoiB3dXJkZSBiZWkgZWluaWdlbiBHZWhhbHRzYW5nYWJlbiDwn5KwIGVudGZlcm50LiAKCmBgYHtyfQppZF9yYW5kb20gPC0gcnVuaWYoODAsIG1pbiA9IDEsIG1heCA9IG4pICU+JSByb3VuZCgpCgpkdCA8LQogIGR0X3JhdyAlPiUgCiAgbXV0YXRlKGVpbmtvbW1lbiA9IGlmX2Vsc2UoaWQgJWluJSBpZF9yYW5kb20sIHJvdW5kKGVpbmtvbW1lbi8xMCksIGVpbmtvbW1lbikpCgpnZ3Bsb3QoZHQsIGFlcyh4ID0gYWx0ZXIsIHkgPSBlaW5rb21tZW4pKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuMikgKwogIHRoZW1lX3N0YXRhKCkKYGBgCgpXaXIgZXJoYWx0ZW4gZWluZSB2ZXJ6ZXJydGUgU2Now6R0enVuZy4g8J+YnwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKcGxvdF9sbShmaXQsIGJ3PTUwKQpgYGAKCkRpZSBmYWxzY2hlbiBCZW9iYWNodHVuZ2VuIG3DvHNzZW4gd2lyIGVudGZlcm5lbiAob2RlciBrb3JyaWdpZXJlbikuIPCfmIMKCmBgYHtyfQpkdF90bXAgPC0gZHQgJT4lIGZpbHRlcihlaW5rb21tZW4gPiA1MDApCgpmaXQgPC0gbG0oZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdF90bXApCnRpZHlfbG0oZml0KQpnbGFuY2VfbG0oZml0KQpwbG90X2xtKGZpdCkKYGBgCgoKCiMjIEF1dG9rb3JyZWxhdGlvbgoKIVtdKGJpbGRlci90ZWxlZm9uemVudHJhbGUuanBnKQoKRGllIFBvc3RiZWFtdGVuIPCfk68gd3VyZGVuIDIwMTAgdW5kIDIwMjAg8J+ThSBuYWNoIGlocmVtIEVpbmtvbW1lbiDwn5KwIGJlZnJhZ3QuIOKYju+4jwoKYGBge3J9CmR0X3RtcCA8LQogIGR0X3JhdyAlPiUgCiAgbXV0YXRlKGVpbmtvbW1lbiA9IHJvdW5kKGVpbmtvbW1lbiAtIHBhcmFtX2FsdGVyKjEwICsgcm5vcm0obiwgc2QgPSAxMCkpLAogICAgICAgICBhbHRlciA9IGFsdGVyIC0gMTAsCiAgICAgICAgIGJlZnJhZ3VuZyA9IDIwMTApICU+JSAKICBmaWx0ZXIoYWx0ZXIgPj0gMjApCgpkdCA8LQogIGR0X3JhdyAlPiUgCiAgbXV0YXRlKGJlZnJhZ3VuZyA9IDIwMjApICU+JSAKICByYmluZChkdF90bXApICU+JSAKICBhcnJhbmdlKGlkLCBiZWZyYWd1bmcpICU+JSAKICByZWxvY2F0ZShiZWZyYWd1bmcsIC5hZnRlciA9IGlkKQoKRFQ6OmRhdGF0YWJsZShkdCAlPiUgZmlsdGVyKGlkID4gNTAwKSkgICMgRGF0ZW4gYWIgSUQgNTAwIHVtIFBhbmVsLVN0cnVrdHVyIHp1IHZlcmRldXRsaWNoZW4gCmBgYAoKUG9zdGJlYW10ZSDwn5OvIGRpZSAyMDEwICDwn5OFIMOkbHRlciBhbHMgNTAgd2FyZW4gc2luZCBpbnp3aXNjaGVuIGluIFBlbnNpb24g8J+RtSB1bmQgaGFiZW4gYW4gZGVyIDIwMjAgQmVmcmFndW5nIOKYju+4jyBuaWNodCBtZWhyIHRlaWxnZW5vbW1lbi4KCmBgYHtyfQpkdF90bXAgPC0gZHQgJT4lIG11dGF0ZShiZWZyYWd1bmcgPSBmYWN0b3IoYmVmcmFndW5nKSkKCmdncGxvdChkdF90bXAsIGFlcyh4ID0gYWx0ZXIsIHkgPSBlaW5rb21tZW4sIGNvbG91ciA9IGJlZnJhZ3VuZykpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC4yKSArCiAgdGhlbWVfc3RhdGEoKQpgYGAKCkltIE1vZGVsbCAobiA9IDE3NTApIHNpbmQgQmVvYmFjaHR1bmdlbiAodW5kIFJlc2lkdWVuKSBuaWNodCBtZWhyIHVuYWJow6RuZ2lnLiDwn5ifIFdpciBoYWJlbiBqYSB2aWVsZSB6d2VpbWFsIGJlZnJhZ3QuIOKYju+4jwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0X3RtcCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCiMgQXVmZ2FiZQoKIVtdKGJpbGRlci9yYWRpb3F1aXouanBnKQoKRGlza3V0aWVyZW4gU2llIPCfmYfwn5mH8J+ZhyBkaWUgUmVncmVzc2lvbnMtQW5uYWhtZW4gIPCfk4ggIHVuZCB6ZWlnZW4gU2llIG9iIHVuZCB3aWUgVmVybGV0enVuZ2VuIPCfj6UgZGVyIEFubmFobWVuIGVya2FubnQg8J+Vte+4jyAgd2VyZGVuIGvDtm5uZW4uCgotLS0KCiMgUXVlbGxlbgoKKyA8aHR0cHM6Ly9jb21tb25zLndpa2ltZWRpYS5vcmcvd2lraS9GaWxlOlNvcGhpZUFuZGVyc29uVGFrZXRoZWZhaXJmYWNlb2ZXb21hbi5qcGc+CisgPGh0dHBzOi8vY29tbW9ucy53aWtpbWVkaWEub3JnL3dpa2kvRmlsZTpCdW5kZXNhcmNoaXZfQmlsZF8xODMtNjMxNDItMDAwMSxfQ290dGJ1cyxfQnJpZWZ0ciVDMyVBNGdlcmlubmVuLmpwZz4KKyA8aHR0cHM6Ly9jb21tb25zLndpa2ltZWRpYS5vcmcvd2lraS9GaWxlOkJ1bmRlc2FyY2hpdl9CaWxkXzE4My0xOTkwLTA5MjQtMDIwLF9HZXJhLF9Qcm90ZXN0X3Zvbl9Qb3N0LUdld2Vya3NjaGFmdGVybi5qcGc+CisgPGh0dHBzOi8vY29tbW9ucy53aWtpbWVkaWEub3JnL3dpa2kvRmlsZTpCdW5kZXNhcmNoaXZfQmlsZF8xODMtVTAyMDUtMDAxMCxfQmVybGluLF9UYWdfZGVzX1Bvc3QtX3VuZF9GZXJubWVsZGV3ZXNlbnMsX1ZvcmJlcmVpdHVuZy5qcGc+CisgPGh0dHBzOi8vY29tbW9ucy53aWtpbWVkaWEub3JnL3dpa2kvRmlsZTpCdW5kZXNhcmNoaXZfQmlsZF8xODMtNTkxMTEtMDAwMixfQmVybGluLF9Qb3N0YW10X05fNTgsX0JyaWVmdHIlQzMlQTRnZXIuanBnPgorIDxodHRwczovL2NvbW1vbnMud2lraW1lZGlhLm9yZy93aWtpL0ZpbGU6VGVsZWZvbnplbnRyYWxlX2JlaV9kZXJfMi5JbmZhbnRlcmllZGl2aXNpb25fKEJpbGRJRF8xNTcwNTEyNSkuanBnPgorIDxodHRwczovL2NvbW1vbnMud2lraW1lZGlhLm9yZy93aWtpL0ZpbGU6UmFkaW9xdWl6XyUyMkFsbGVpbl9nZWdlbl9hbGxlJTIyX2F1c19kZW1fUmF0c3NhYWxfKEtpZWxfNDUuNTY3KS5qcGc+CgotLS0KCiFbXShiaWxkZXIvZmFpci5qcGcp